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Abstract 

A new numerical treatment in the Crank-Nicholson method with the imaginary time evolution 
operator is presented in order to solve the Schrodinger equation. The original time evolution 
technique is extended to a new operator that provides a systematic way to calculate not only 
eigenvalues of ground state but also of excited states. This new method systematically produces 
eigenvalues with accuracies of eleven digits with the Cornell potential that covers non-perturbative 
regime. An absolute error estimation technique based on a power counting rule is implemented. 
This method is examined with exactly solvable problems and produces the numerical accuracy 
down to 10~^^. 



1 



Numerical computation of analytically unsolvable Schrodinger equations has been of in- 
terest in atomic and molecular physics, quantum chromodynamics, and Bose-Einstein con- 
densation of trapped atomic vapors . Conventionally, a wave function has been 
represented as a linear combination of plane waves or of atomic orbitals [5]. However, these 
representations entail high computational cost to calculate the matrix elements for these 
bases. The plane wave bases set is not suitable for localized orbitals, and the atomic orbital 
bases set is not suitable for spreading waves. In particular, potential problems such as the 
Cornell potential 0, ^ are difficult to compute precise eigenvalues because they have sin- 
gularities at the origin and at the infinity, and includes non-perturbative regime when the 
linear term is significant. 

To ove.c„.e these „„b.e,„s. ,.,„e*a. .ethod. such as Land. sub.ae«o„ and Nys- 
trom plus correction ^ in the momentum space produced eigenvalues with six or seven 
digits. Others adopted real-space representation 10,[ll|. In these methods, a wave func- 
tion is discretized by grid points in real space providing from five to seven digits. Also, 
estimates using the exact solutions of the Killingbeck potential produced eigenvalues with 
seven digits Among these real-space methods, a method called Crank-Nicholson (C-N) 
scheme 12|, |l3| is known to be especially useful for one-dimensional systems because this 
method conserves the norm of the wave function exactly and the computation is stable and 
accurate even in a long time slice. These characteristics are very attractive for solving the 
Cornell potential in order to compute eigenvalues and eigenstates of the system precisely. 

The current numerical precisions in relevant research areas may not be as high as one 
would hope to achieve. For example, it may forbid one to study fine or hyper-fine structure 
effects in the atomic system Calculation of matrix elements subiect to large subtraction 
may require high accuracy in quantum chromodynamics jlj,|15|. Also, none of the references 
we compiled for this study contains the serious error estimate on their numerical calculations, 
which is an important indicator of the reliability of a suggested numerical method. 

In this letter, we apply the C-N method to solve the Schrodinger equation and present 
two new numerical methods. First, the C-N method with the imaginary time evolution 
operator is re-interpreted by extending its allowed region. This method produces ground- 
state eigenvalues with numerical accuracies of eleven digits when the Cornell potential is 
used. We then extend the original time evolution technique to a new operator that provides 
a systematic way to calculate not only the eigenvalue of the ground state but also those 
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of excited states with less computational time by a factor of 10, while providing the same 
accuracy as in our re-interpreted C-N method. Our methods will be useful in calculations 
of eigenstates for the Cornell potential and we discuss some of the results. At the end, we 
discuss a mathematically simple but rigorous absolute error estimation on the numerical 
calculations presented in this letter. 



C-N method 



12l llS^ is a finite difference method used for solving diffusive initial- value 



problems numerically. The main part of this method is the time evolution operation and t 
evolution operator for the Schrodinger equation may be approximated by Cayley's form 
as 

-mt _ ^ ~ 2'^^ I mfn-/3+3 



he 

B 



= Y:d^ + o(HV), (1) 

where Ti. is the Hamiltonian of the problem. This equality is correct up to the second order 
in Tit and the approximation is valid when |7it| ^ 1. By multiplying this operator to an 
initial wave function, one obtains the evolved wave function. The standard C-N method 
makes use of Eq. ((T)) in order to study the time evolution of the wave function jl^ l]^. 



Next, we introduce the imaginary time method 



a 



to calculate the eigenfunctions and 



eigenvalues. By a Wick rotation, one replaces t by —ir in the time evolution operator j^. 
This transforms the original Schrodinger equation into a diffusion equation. Then a wave 

function evolves in time slice At as 

00 00 
u{p, At) = ^ CMp)e-'^''^' = Yl CMp)e-^''^\ (2) 

1=1 i=l 

where Ui{p) and Q are the eigenfunction and the eigenvalue for z-th state, respectively. Q 
is the relative amplitude for i-th state and the summation is over all possible eigenstates of 
the system. 

For the imaginary-time version, eigenfunctions decay monotonically in time till the steady 
state is reached. The ground-state eigenvalue can then be read off from the steady-state 
eigenfunction as r ^ 00 0|. Therefore, we acknowledge that the time evolution operation 
itself in the C-N method acts as a tool that selects the ground state only. Here, we advocate 
that the condition |7Yr| ^ 1 is not necessary in the calculation of the ground-state eigenvalue. 
When all the eigenvalues are negative such as the pure Coulomb potential case, where 
Co < ■ ■ ■ < Cn < • ■ • < 0, the amplification of the ground-state coefficient may happen in 
the region —2 < Tir < as the time evolution continues. On the other hand, when all the 
eigenvalues are positive such as in problems with the linear and Cornell potentials, where 



< Co < ■ ■ ■ < Cn, the time evolution operator can amplify the ground-state coefficient in 
the region Tir < —2. Note that we relax the condition |7Yr| ^ 1 in order to obtain faster 
convergence to ground-state eigenfunction even if that region is used in the standard C-N 
method. We call this approach as a relaxed C-N method. 

Once the ground-state wave function is obtained, one can obtain the ground-state eigen- 
value from the expectation value of the Hamiltonian of interest. Note that, for the numerical 
computation of an expectation value, the upper bound of the integral can not be an infinity 
but a cut-off value, Pmax- We will explain how to control the numerical error produced by 
ignoring the region (pmax, oo) later. 

The spherical Stark effect in hydrogen and a bound state for a heavy quarkonium may 
be described by the Schrodinger equation with the Cornell potential Q, Bl • The radial 
part of the Schrodinger equation with orbital-angular-momentum quantum number i that 
is relevant to such problems may be given by 

£{£ + 1) A 



u{p) = Cu{p), (3) 



dp"^ p2 p 

where the dimensionless wave function u{p) and the dimensionless energy eigenvalue ( in 
Eq. (jS)) are described in Ref. The parameter A is the relative strength between the 
Coulomb and the linear potentials. We call it as the coulombic parameter throughout this 
letter. 

The relaxed C-N method is now applied to the Cornell potential problem. The A de- 
pendence of the ground-state eigenvalues in Ref. jd| is reproduced with our relaxed C-N 
method and a comparison of two results is summarized in Table H] The time evolution in 
our relaxed C-N method gives eigenvalues over iterations giving stable sixteen-digits and 
all agree well with the results in Ref. jfl]. We also find that the convergence speed is im- 
proved by 10 times compared with the standard C-N method that we tested. For excited 
states, one can in principle obtain the eigenvalues from the lowest to higher states by the 
Gram-Schmidt orthogonalization procedure. We find that the standard C-N method is not 
useful in calculating excited states because the convergence speed is practically zero when 
we require the numerical precision to be 10 digits. However, with the relaxed C-N method, 
the excited-state eigenvalues are successfully found, with the speed being 10 times slower 
than the time needed to find the ground state. The amount of time required to calculate 
excited states is no longer than ten minutes with a Pentium IV CPU with 0.5 GB memory. 
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TABLE I: Dependence of the ground-state eigenvalues C on the coulombic parameter A. The 
Cornell potential is used in the calculation. The second column contains the numerical results 

n 

from Ref. [o| and the third column contains our result with the relaxed C-N method. The number 
of grid points is set to be = 300,000 and Ar = — 5 for the numerical analysis. An estimate on 
the numerical errors of the computation AC is also listed at the last column and will be discussed 
later. 

We emphasize that the condition |7^r| ^ 1 can be relaxed as far as the goal of the 
numerical computation is to obtain the ground-state eigenfunction. We extend this idea 
further in order to calculate excited-state eigenvalues and eigenfunctions systematically and 
more efficiently. 

Consider the following operator 

where (3 is an arbitrary real number in the eigenvalue space. If we apply the operator in 
Eq. (jH) to the wave packet k times, 

k ^ 

:Ui{p)^ (5) 



(0 - P)' 



where we assume that the time independent wave packet u{p) can be expressed as a linear 
combination of the infinite number of eigenfunctions Ui{p) with coefficients Ci. For (3 = Q + e 
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such that |e| ^ 1, we have 



1 




1 




> 




0-/5 


0-/5 



(6) 



if j 7^ i. Therefore, Eq. Q plays a role as an amplifier which magnifies the contributions of 
the term with the nearest eigenvalue from the point (3. In this way, all the eigenvalues within 
an arbitrary range in /3 can be found systematically, by running (3 within the range. We call 
this approach as a modified C-N method. Advantages of this method are as follows. First, the 
calculations of excited states can be carried out systematically by stepping through different 
values of (3. Second, the computing time for calculating excited-state wave functions or 
eigenvalues is similar to that needed for the ground state. Furthermore, it does not lose 
accuracies in the calculation of higher-state eigenvalues while other methods do often. This 
contrasts the modified C-N method to the relaxed C-N method. With the relaxed C-N 
method, the Gram-Schumidt orthonormalization slows down the computing speed by ten 
times, as mentioned before. 

First, the modified C-N method is tested on the pure Coulomb potential and to the 
pure linear potential (A = 0) where the exact eigenvalues are known for both cases. This 
is a good benchmark because one directly examines the performance of the algorithm by 
comparing the exact solutions to the numerical results of the algorithm to be tested. We 
find that the numerical values of eigenvalues agree well with known analytical values, up 
to eleven digits. Second, the modified C-N method is applied to the Cornell potential and 
reproduce the result in Table HI We find that eleven digits of the eigenvalues are reproduced 
completely from A = 0.0 to 1.8. This is consistent with our error estimation which will be 
explained later. Third, the modified C-N method is applied to the calculation of excited 
states. Table |n] shows the eigenvalues obtained. Note that for the IS* state in Table m 
one can compare the eigenvalue with the number in Table H] and only the last two digits 
are different, which is again consistent with our error estimation. Again, for the ground 
state, the relaxed and the modified C-N methods require similar amount of computing time, 
but for the excited states, the modified C-N method is faster by 10 times. We checked the 
computing speed for the Coulomb, linear, and Cornell potentials separately and all three 
give similar performances. 

There are two sources of errors in our numerical calculation. One is the cut-off (pmax) 
and the other is the discretization of continuous equations. The cut-off gives the imperfect 
numerical integration but could in principle be reduced as small as one wishes by increasing 
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TABLE II: Numerical values of eigenvalues and the error estimation for various excited states with 
the Cornell potential. The modified C-N method is used to calculate eigenvalues. The coulombic 
parameter is set to A = 1.0 and the number of grid points to be AT = 300, 000. 

the value of pmax- We estimate the error due to the finite value of pmax for the Coulomb 
and Unear potentials, respectively, by integrating exact eigenfunctions. Wc find that the 
error is 10~^^ or smaller when Pmax = 20, for example. In fact, for the practical purpose, we 
control the value of pmax in such a way that the numerical error due to finite value of pmax is 
smaller than the error from the fact that we deal with discrete processes. This assumes that 
the error estimates with the Coulomb and linear potentials individually are not significantly 
different from the errors due to the Cornell potential. Practically, selecting proper values 
of Pmax is important, because too small values cause large errors and too large values may 
slow down the computation. Wc find that Pmax = 30 is an optimal value for most of our 
applications shown in this letter. Note that pmax — 30 corresponds to 30 times of the Bohr 
radius in the hydrogen atom problem, for example. 

More serious source of the error is originated from the discretization of continuous differen- 
tial equations. In this letter, differentiation and integration of wave functions are discretized 
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FIG. 1: Numerical error estimation for the eigenvalues as a function of the number of grid points. 
For the Coulomb (linear) potential, ICexact — C(-^)l is indicated as open boxes (triangles) for IS, 
IP, and ID states (IS* state only for the linear potential). The error estimates based on Eq. © for 
the Coulomb, linear, and Cornell potentials are shown in straight, dashed, and dot-dashed lines, 
respectively. Note that both axes are in logarithmic scales. 

with the following prescriptions 

/ dpuip) = ("^+1 + + ^( V), (7b) 

where Ap is the distance between two nearest discrete points in the calculation. From 
both Eqs. (fTaj) and ()7b|) . numerical errors contained in the discretization are proportional 
to IS.f? = A^~^ where N is the number of grid points in the discretization. Therefore, the 
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error in eigenvalues due to the discretization may be approximated as 

AC(iV) ^ (8) 

where AC(A^) = iCexact — C(^)l the constant c depends on the potential. If we select 
two arbitrary values in grid points, A^i and for example, then we can easily obtain the 
constant c as 



With this, one can estimate the error due to the discretization prescription in Eqs. (fTaj) 
and ()7b|) . We refer it as an error estimation from the power counting rule. In Fig. Q our 
estimate of the numerical error A<^(A^) is compared with the true error ICexact — C(^)l for the 
Coulomb and linear potentials. Here we used our modified C-N method for the error analysis. 
It is apparent that our error estimate is accurate down to 10~^^ for the IS" state under the 
Coulomb potential, for example, when iV=300,000. For others, the results are better as in 
Fig.^ For the ID state under the Coulomb potential, the true error looks unstable at large 
value of and we find that this is due to the limitation in storing significant digits during 
our computation of C{N). Note that for the Cornell potential, the true errors cannot be 
calculated because the exact solutions are unknown. Therefore, for the Cornell potential, we 
estimate that the errors are in the range of 10~^^ as in Fig. ^ We use this error estimation 
technique throughout this letter and numerical values of the error estimates are included in 
Tables m and nil 

In conclusion, we have presented two numerical methods for calculating Schrodinger 
equation with the Crank-Nicholson method. In the relaxed C-N method, the time evolution 
operator is re-interpreted as a weighting operator for finding the ground state eigenfunction 
more precisely. This idea is extended to a new operator in the modified C-N method that is 
more efficient in computing not only the ground-state but also excited-state wave functions 
systematically. An absolute error estimation method is presented based on a power counting 
rule and is consistent with predictions when exact solutions are known. These two algorithms 
may be useful when precise numerical results are required. Possible applications may include 
Cornell potential 0,Q] and Bose-Einstein condensation of trapped atomic vapors ^. 
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